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ABSTRACT 

The source of viscosity in astrophysical accretion flows is still a hotly debated issue. We 
investigate the contribution of convective turbulence to the total viscosity in a self-consistent 
approach, where the strength of convection is determined from the vertical disc structure itself. 
Additional sources of viscosity are parametrized by a /3-viscosity prescription, which also 
allows an investigation of self-gravitating effects. In the context of accretion discs around 
stellar mass and intermediate mass black holes, we conclude that convection alone cannot 
account for the total viscosity in the disc, but significantly adds to it. For accretion rates up 
to 10% of the Eddington rate, we find that differential rotation provides a sufficiently large 
underlying viscosity. For higher accretion rates, further support is needed in the inner disc 
region, which can be provided by an MRI-induced viscosity. We briefly discuss the interplay 
of MRI, convection and differential rotation. We conduct a detailed parameter study of the 
effects of central masses and accretion rates on the disc models and find that the threshold 
value of the supporting viscosity is determined mostly by the Eddington ratio with only little 
influence from the central black hole mass. 
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. 1 INTRODUCTION 

Modern theoretical modeling of accretion discs dates back to the 
year 1948, when Weizsacker published his article about the rotation 
, of cosmic gas (Weizsacker 1948). A key ingredient to describe the 
' accretion process is the origin of the viscosity, which causes fric- 
tion in the disc and an inward motion of the material. First obser- 
vations of accretion timescales of discs in cataclysmic variables in- 
validated theoretical expectations of molecular viscosity being the 
, driving force (see, e. g., Prendergast & Burbidge 1968; Pringle & 
< Rees 1972). Adversely, they revealed a discrepancy of many orders 
of magnitude between the numbers measured in the lab and those 
needed to account for the observations. Soon thereafter, Shakura 
& Sunyaev (1973) proposed the a-viscosity parametrization, by 
which most observations could be reproduced satisfactorily. Never- 
theless, the Shakura-Sunyaev viscosity remains a purely empirical 
description and is limited to thin discs with negligible disc masses 
(i.e., non-selfgravitating discs). Among the physical theories, the 
most promising ones are: 

Differential rotation. An obvious candidate for the turbulence 
in nearly Keplerian rotating discs is differential rotation. From early 
laboratory experiments on rotating Couette-Taylor flows (Wendt 
1933; Taylor 1936), this possibility was first ruled out. However, 
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in recent re-investigations, Richard & Zahn (1999) and Richard 
(2001) concluded that differential rotation can give rise to tur- 
bulence, despite published arguments. At the same time, Duschl, 
Strittmatter & Biermann (1998, 2000) formulated the viscosity 
description. Although being a parametrization like its ancestor, it 
can actually be related to the process of differential rotation. Con- 
trary to the Q-prescription, the /3-viscosity accounts properly for the 
selfgravity of the disc. At the same time, it includes the a-viscosity 
in the case of a shock dissipation limited, non-selfgravitating disc. 
Combining the laboratory measurements with the formulation of 
the /3-viscosity leads to a value of /? ~ 10~^ which can be pro- 
vided by differential rotation.^ 

Convection. In order to account for the transport of the en- 
ergy released by the accretion process, convection is considered to 
support or even dominate in some cases over radiation and has been 
studied intensively (Bisnovatyi-Kogan & Blinnikov 1977; Shakura, 
Sunyaev & Zilitinkevich 1978; Goldman & Wandel 1995; Agol et 
al. 2001) with substantially different conclusions (not least due to 
the underlying theoretical models): the contribution of convection 
to the overall energy transport regions ranges from 1/3 (Shakura 
et al. 1978) to being completely dominant (Bisnovatyi-Kogan & 
Blinnikov 1977) in radation pressure dominated disc regions. In 

^ As a rale of thumb, (3 ^ . . .a (Duschl et al. 2000); the corresponding 
a parameter therefore lies between 10~^ and 10~^. 
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gas pressure dominated regions, multiple solutions are found in the 
same range (Goldman & Wandel 1995). Also, recent 2-dimensional 
simulations by Agol et al. (2001) demonstrate that convective pro- 
cesses can release heat sufficiently fast to modify the vertical struc- 
ture of the disc. It is therefore natural to consider the turbulence 
caused by convective motion as a possible candidate for viscosity. 
First (semi-)analytical investigations were discouraging: they led 
to discs with masses comparable to or even exceeding the central 
black hole mass (Vila 1981; Duschl 1989), incompatible with the 
a-viscosity description assumed in the models. Ruden et al. (1988) 
and Ryu & Goodman (1992) studied convective instabilities in thin 
gaseous discs and confirmed that angular momentum transport can 
be supported by convective turbulence. Goldman & Wandel (1995) 
investigated accretion discs where viscosity is given by convection 
solely and where the energy transport is maintained by radiation 
and convection. They found the resulting viscosity being too low 
by a factor of 10 to 100, but could not draw final conclusions due 
to their limited disc model. 

Magneto-rotational instability. The magneto-rotational insta- 
bility (MRI) was first noticed in a non-astrophysical context by 
Velikhov (1959) and Chandrasekhar (1960). More than 30 years 
later, Balbus & Hawley (1991, 1998) established that weak mag- 
netic fields can substantially alter the stability character of accre- 
tion discs, giving rise to a generic and efficient angular momentum 
transport. Today, the MRI is considered as the primary candidate 
for the viscosity in astrophysical accretion flows. Modern compu- 
tational facilities allow the study of angular momentum transport 
in magnetized discs in 3-dimensional MHD codes, which basically 
can be used to calibrate the a- or /3-viscosity parameter (see, e. g., 
Balbus 2005, for a review). The key problem therein is the non- 
trivial dependency of a on various physical and numerical param- 
eters of the simulations. Recently, Pessah, Chan & Psaltis (2007) 
presented a scaling law which allows to disentangle physical and 
numerical influences. The general question if the MRI effects can 
be translated into an a- or l3-\ypt viscosity remains to be answered 
(see, e. g., Pessah, Chan & Psaltis 2008). From the wealth of results 
obtained so far, it seems likely that the MRI alone cannot account 
for the viscosity in astrophysical discs (Begelman & Pringle 2007; 
King et al. 2007). For example, current results face a discrepancy of 
at least one order of magnitude between the viscosities generated 
by the MRI and those inferred from observations (see also Lcsur 
& Longaretti 2007, for a further discussion). A particular problem 
of the MRI are the so-called dead zones, where the growth rate of 
magneto-rotational instabilities is strongly suppressed and the tur- 
bulence induced by magnetic effects diminishes (Gammie 1996). 
Although the implications of MRI dead zones are discussed mostly 
for protoplanetary discs (see, e. g., Reyes-Ruiz, Perez-Tijerina & 
Sanchez-Salcedo 2003; Brandenburg 2(X)8, for an overview), the 
overall problem of a vanishing viscosity applies to accretion discs 
in general. 

Thus, it is not yet clear whether one of these candidates or 
a combination of them is responsible for generating the viscos- 
ity in astrophysical discs. One important step therefore is to study 
the effect of convective turbulence in combination with other con- 
tributors. Goldman & Wandel (1995) stressed the need for a con- 
vective disc model where the vertical structure is calculated self- 
consistently in order to quantify better the convective turbulence 
and the energy transport in the disc. 

In this paper, we construct a model of a black hole accretion 
disc where we calculate the effect of convection in a self-consistent 
way by means of the mixing-length theory. Hereby, the total viscos- 
ity is given by convection plus a supporting /?- viscosity, accounting 



for turbulence due to differential rotation and allowing for potential 
self-gravitating effects. Energy transport in the vertical direction is 
provided by radiation and convection simultaneously, which allows 
to derive the strength of the convective viscosity within the model. 
The details of the model are given in Sect. 2; in Sect. 3, we present 
and analyse our results for various central black hole masses, ac- 
cretion rates and values of the underlying /9-viscosity. Section 4 is 
devoted to discussion and conclusion. 



2 MODEL SETUP 

We calculate accretion disc models where both the viscosity and the 
transport of energy is supported by convective processes, in addi- 
tion to an underlying /3-viscosity and to radiative energy transport. 
Hereby, convection is treated in the framework of the mixing-length 
theory. The disc is assumed to be geometrically thin in order to al- 
low for a 1 -|- 1-dimensional treatment of the equations. We use a 
cylindrical coordinate system with planar radial coordinate s, verti- 
cal coordinate z and true radius r = Vs^ + 2^- The disc geometry 
is determined by an inner and an outer radius, Si and So. and the 
disc's thickness h = h{s) from the mid-plane. 

The turbulent viscosity, caused by convective processes, is 
generally given by 

fcom ~ ^^conv^conv (1) 

with Vconv and /com being the turbulent velocity of the convective 
elements and the convective lengthscale over which they diffuse, 
respectively. The factor ^ is of the order of unity and depends on 
the degree of isotropy of convection in the considered direction. For 
simplicity, we assume isotropy in this investigation (i. e., ^ = 1/3). 
We identify the convective lengthscale Zconv with the mixing-length 
lia, which will be deflned later. We include other sources of vis- 
cosity (differential rotation, MRI, . . . ) by assuming a permanently 
supporting viscosity to be present in the disc, parametrized by a 
standard /3-ansatz: 

= fis'uj , /3 < 1 . (2) 

Here, oj stands for the angular velocity. The total viscosity v is then 
given by a combination of these two contributors, 

V = Vcom + Vf} ■ (3) 

2.1 Radial structure 

For the calculation of the radial structure, we introduce 

* = A vpdz. (4) 
Jo 

in analogy to the disc's surface density S, 

E= f pdz. (5) 
Jo 

Only in the special case ofv = const, can we rewrite (4) to = 
vT,. In all other cases, we apply the mean value theorem to define 
an average value v* such that 

* = i>*T, 

The individual contributors v* and E remain unknown from the 
radial structure equations only, but are determined by the vertical 
structure equations (Sect. 2.2). 
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The radial structure is determined by the conservation of mass, 
momentum, angular momentum and energy. The corresponding 
equilibrium equations are 



OJ 

2* 
2F 



s 



Mui 



2ns{du}/ds) 

Msoj{duj/ds) 
27r 




M denotes the (constant) accretion rate, Vs the accretion velocity 
with Vs <Q for inflowing material, Qs the graviational acceleration 
in radial direction, and F the heat flux, integrated in vertical direc- 
tion. We apply the standard free-fall boundary condition (Shakura 
& Sunyaev 1973; Novikov & Thorne 1973) in (8) and (9), imply- 
ing a vanishing torque at the disc's inner radius s\. The momen- 
tum equation (7) is simplified by the assumption of local equilib- 
rium of the graviational attraction and the centrifugal repulsion in 
the radial direction. Therein, the gravitational acceleration is as- 
sumed to be given by the monopole approximation (Mineshige & 
Umemura 1997), assuming a Pseudo-Newtonian gravitational po- 
tential (Paczyiiski & Wiita 1980): 



G{M^ + Mi{s)) _ s 
{r — rs)^ r 



(10) 



The Schwarzschild radius rg is given by 2GMc/c^, where Mc de- 
notes the mass of the central black hole. The enclosed disc mass at 
radius s is calculated via 



-f 

J Si 



47rs'S ds . 



(11) 



2.2 Vertical stratification 

2.2.1 Structure equations 

In analogy to Cannizzo & Cameron (1988), Hofmann (2005) and 
Vehoff (2005), we adopt the energy flux at height z as the indepen- 
dent coordinate for the vertical integration: 

-l^dz, F,{z = h)=F. (12) 
Jo oz 

Additionally, we introduce the surface density at height z, 

E, = / pdz, E,(2 = /i) = S, (13) 

^0 



and 



■0 = / vpdz, = = 

Jo 



(14) 



Neither E nor h are known a priori - they will be a result of the 
vertical integration. The equations for the vertical structure of the 
disc are given as follows: 



dz 



(6) 


or 




(7) 




(8) 


9E^ 




dF, 


(9) 


dp 







1 



(15) 

(16) 
(17) 

(18) 



Here, Cp stands for the isobaric specific heat capacity and 7 = 
Pgm/p- (15) is a simple inversion of the local energy production by 
viscous dissipation, dF^/dz = pi's^{dLo/ds)^ . The temperature 
stratification (16) results from accounting for the energy transport 
by radiation and convection (Cox & Giuli 1968) and using (15). Re- 
lating the two terms in (16) with the radiative and adiabatic gradi- 
ents Vrad and Vad, it becomes clear that the variable describes the 
relative contribution of the convective energy transport to the total 
energy transport. Its value depends on the local physical conditions 
at position (s, z) in the disc and can be calculated numerically, see 
Sect. 2.2.2. Combining (15) with the definition of the surface den- 
sity at height z (13) leads to (17). Assuming hydrostatic equilibrium 
and again using (15) gives the last differential equation (18) for the 
pressure stratification. 

Using (16), (18), and the equation of state. 



pkjiT 4(TSB„4 
P = Pgas + Prad = + , 

we transform (18) into an equation for the mass density p: 

4asBr" 



(19) 



dp_ 

OF, 



pc 



4acT3 



7Cp 



(20) 



For the numerical solution of the vertical stratification, the opacity 
K = k{p, T) is calculated from a combination of tabulated values 
and interpolation formulae, see Sect. 2.2.4. In analogy to Qs, the 
gravitational acceleration in vertical direction Qz is provided by the 
monopole approximation: 



G(Mc + M,{s)) z 

9z = 7 • - + 47rGEj . 

(r — rs)^ r 



(21) 



In (21), the second term stands for the local gravitational attraction, 
which becomes important in the self-gravitating regime. 



2.2.2 Adaptation of the mixing length theory 

In order to solve the vertical structure equations, we apply the 
mixing-length theory (Bohm-Vitense 1958) as formulated in Cox 
& Giuli (1968). The mixing-length theory expresses the efficiency 
of the convective energy transport relative to the radiative transport 
processes by the variable C, where ^ ^ ^ 1. A vanishing ^ im- 
pUes no convective transport, while in the case = 1 all energy is 
transported by convection. Following Cox & Giuli (1968), its value 
can be calculated from the cubic equation 



^1/3 ^ g . ^2/3 _^ ^^^2^ _ ^^^2 ^ Q ^ 



(22) 
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with a numerical factor ao = 9/4 and further quantities defined as 

1 1/3 

Vad) 



(Vrad • 



Q 



Cp = 



288a2c2pT6 
4-37 



7 



32 - 247 - 37^ 



272 



Vrad 
Vad 



4acT4 ' 
8-67 



32 - 247 - 372 ' 
min(Ap, h) , 

Ap — . 

The mixing-length Im is usually of the order of the pressure scale 
height Ap. However, in analogy to the stellar case, it is limited by 
simple geometric effects. While in the stellar case, it usually can- 
not exceed the actual radial distance from the centre due to sym- 
metry requirements, we adopt the actual height h = h{s) of the 
disc as upper Umit. In doing so, convective elements are allowed to 
travel across the disc mid plane, which overrides the symmetry of 
the disc. Also, it removes the strict upper barrier (i. e., the disc sur- 
face) for the convective elements. Nevertheless, it provides a simple 
method of taking into account overshooting effects and a more re- 
alistic, smooth transition between the disc and the atmosphere. We 
note that the results differ only slightly for a more restrictive defini- 
tion = min(Ap, h — z, z), so that our conclusions do not depend 
on this assumption. 

From the above definitions, the cubic equation (22) is solved 
numerically. Subsequently, the convective viscosity is calculated 
from (Cox & Giuh 1968) 



aoA 



1/3 



c 



1/3 



(23) 



where Cs denotes the sound speed at the actual coordinate (s, z) in 
the disc. In the non-relativistic regime, it is given by 



(24) 



The constant Fi stands for the poly tropic index, which is given by 
5/3 in the case of a non-relativistic, ideal gas. The parameter q;„ 
relates the typical distance traveled by the convective elements to 
the pressure scale height: 



An 



(25) 



The derivation of the mixing-length theory assumes a purely sub- 
sonic motion of the convective elements, Vconv < Cs. This inequal- 
ity cannot be assured by the definition of the convective velocity 
in (23). Hence, in the case (23) leads to values Vcom > Cs, we 
follow Cox & Giuli (1968) and set manually iiconv = Cs. This, in 
turn, invalidates equation (22) for ^. Instead, ^ must be calculated 
from (23) with 



C = 



-A 



(26) 



Q3/2af J Vrad - Vad) ■ 

This correction leads to lower values of C, and therefore to an effec 
tive decrease of the convective energy transport and of the convec- 
tive viscosity Vcom- 



2.2.3 Boundary conditions 

We define the boundary conditions for the four equations (15)-(17), 
(20) at either the disc's mid plane ("mp", = 0) or surface ("eff", 
= F): 



«mp = z{F:, = 0) = , 
Teff = T{F, =F) = 

VO-SB, 



S^,mp = S^F^ =0) = 0, 



Pe£f = P{Fz = F) = Peff,mpm ■ 



(27) 
(28) 
(29) 
(30) 



The boundary condition on p cannot be determined from the radial 

structure equations or from simple geometric arguments. It rather 
requires us to define an atmosphere above the disc, which allows us 
to determine the density peff.inpui consistently with the height h, the 
surface density E and the effective temperature Teff. Details about 
the atmosphere will be given in Sect. 2.2.5. 



2.2.4 Opacity k 

The opacity has a rather strong influence on the computation of the 
vertical structure. Therefore, we refrain from simple models such 
as pure electron scattering or Kramer's law. We rather use a com- 
bination of tabulated values and interpolation formulae to calculate 
the opacity for a broad range of temperatures and densities. The 
presence of a disc and an atmosphere implies the calculation of 
Rosseland and Planck opacities. 

Tabulated values Given that we want to cover a large domain 
in temperature and density, multiple sources are included in our 
model. In the high-temperature limit, we adopt the tables from the 
TOPS project (TOPS 2008). We compile tables for Rosseland and 
Planck opacities in the range of 

logp = [-12.5 . . . -M0.5] , 
logT = [-F4.5... -1-9.1] . 

All values are given in cgs-units. The number of data points is 47 
on an equidistant scale for log p and 41 for log T, respectively. 

In the low-temperature regime, we include the Ferguson opac- 
ities (Ferguson et al. 2005; Ferguson 2008). We compile Rosseland 
and Planck opacity tables in the range of 

logi^ = [-8.0 . . . -I- 1.0] , 
logT = [-1-2.7... -1-4.5] , 

with R = p/Ti (Te = T/IO"^). These ranges correspond to min- 
imum and maximum mass densities of logp = —17.9 and +6.5, 
with a resolution of 19 equidistant points in log R and 85 in log T. 

We choose identical chemical abundances for the TOPS and 
Ferguson opacities with mass fractions X = 0.7,Y = 0.28, Z = 
0.02, and the chemical mixture of Grevesse & Sauval (1998). 

Analytic interpolation formula Opacities outside the ranges 
given above are calculated using an analytic interpolation formula 
(Gail, priv. comm.; for a similar approach, see also Bell & Lin 
1994) for Rosseland opacities. Thus, strictly speaking, this inter- 
polation formula is vaUd only in the optically thick regions. 
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Table 1. Interpolation of the opacity: set of parameters (in cgs-units). 



Contributor 1 


Symbol 


Kl,0 


Kl,p 


Kl,T 


Dust with ice mantles 




2.0 


10-4 





2 


Evaporation of ice 


'^ice, evap 


1.0 


•IQie 





-7 


Dust particles 


'^dust 


1.0 


10-1 





1/2 


Evaporation of dust particles 


^dust, evap 


2.0 


■1081 


1 


-24 


Molecules 


^mol 


1.0 


10-8 


2/3 


3 


Negative hydrogen ion 


fCjj- 


1.0 ■ 


10-36 


1/3 


10 


Bound-free, free-free-transitions 




1.5 


.1020 


1 


-5/2 


Electron scattering 







348 









T'lO _l_ T'lO 1^4 

J.Q T J- ft-ice, evap 



1/4 



1 



+ '*mol + '*H- 



+ 



^dust, evap 

The individual contributors ki are approximated by 



1/4 



Kl = Ko.l 



P 



(31) 



(32) 



and are compiled in Table 1. The temperature To parametrizes 
the transition between atomic/molecular ice and gas and is set to 
To = 3000 K. Note that the definitions of To and of the individ- 
ual contributors are such that the full interpolation formula (31) fits 
the values obtained from experiments and numerical calculations; 
they cannot be used on their own as a physical descriptions of the 
corresponding processes. 

Opacity mixture For smooth transitions and a broad coverage 
in the T-p range, we use a combination of the tabulated opacities 
(TOPS, Ferguson) and the interpolation formula (31). The transi- 
tion between the TOPS and the Ferguson opacities takes place at 
log T = 4.5, modeled by a linear interpolation of the opacities 
from both sources in the range log T = [4.0 . . . 5.0]. 

At the "outer" boundaries of the TOPS- and Ferguson- 
opacities, we use the same kind of Unear transition in a range of 
A log p = 1 and A log T = 1 between the tabulated values and 
the interpolation formula (31). The resulting opacities are defined 
on a log T-log p grid with 150 data points in each direction and 

logp = [-15.0... 10.0] , 
logT = [-hl.O . . . -I- 9.0] , 

which is sufficient for our purposes. Figures la,b display the result- 
ing opacities as a function of temperature for certain densities. 

2.2.5 Atmosphere 

The only purpose of the atmosphere is to provide a value for the 
mass density at the surface of the accretion disc at each radial po- 
sition s, which is consistent with the actual effective temperature, 
geometrical height and surface density. Therefore, it is sufficient to 
calculate a simple grey atmosphere in the Milne-Eddington way, 
where the temperature distribution is given as a function of the op- 
tical depth r by 



T{t)* = Ci,atmTeff • (t -|- Cg.atm) . 



(33) 



The constants Ci,atm and (72,atm depend on the transition point 
Tetf between the atmosphere (optically thin) and the disc (opti- 
cally thick) and the final value for the temperature at the "upper" 



end of the atmosphere (t <C 1). We use the common value of 
T*{t = 0) = (1/2) • Teff, but allow the transition to take place 
at Teff — 1 (instead of the common value 2/3) for a simple reason: 
in the optical thin atmosphere, the equation of state (19) is modi- 
fied such that the radiation pressure term tends to zero for r ^ 0. 
With the approximate expression for the radiation pressure in an 
optical thin medium (see, e. g., Artemova et al. 1996, for a detailed 
discussion), the atmospheric equation of state becomes 

pfesT 4CTSB „4 
P = Pgas + Prad,atm = h -5 tT . (34) 

Hmn 3c 

A smooth transition of the pressure between the disc and the at- 
mosphere requires Teff = 1, which implies Ci,atm = 1/2 and 
C2,atm = 1. The remaining equations are given as follows: from 
the definition of the optical depth, 

dr — —Kpdz , 

we get an expression for dz/dr. The differential expression for the 
surface density 

dTiz = pdz 

is transformed into 9Ez/9r. We choose the gas pressure as the 
fourth dependent variable and assume hydrostatic equilibrium, 
dp/dz = —pQz, to obtain the following set of differential equa- 
tions for the structure of the atmosphere: 

dz^ 1_ _ fcsT 

dr pK /irriH/tPgas ' 

as, _ _i 

dr K ' 



dT 
dr 



Teff. , - 



9/4 



-3/4 



(35) 
(36) 
(37) 



4(TSB 4 
—5 Jeff • 

3c 



1 

^+2 



(38) 



Equation (37) is derived from the temperature profile (33), and (38) 
is calculated from 



apgas _ dp C)prad,atm _ dz 4asB d 



(rT*) . 



dr dr dr 
The corresponding boundary conditions need to be set at either 
the lower boundary (i. e., at the disc surface, corresponding to r = 
Teff = 1) or the upper boundary ("up", corresponding to t = Tup 
1). Three of these boundary conditions are provided by the solution 
of the vertical disc structure: 



h = z{Fz =F) = z{t = Teff) , 

E = = F) = E,(T = Trff), 

Teff = T{F, =F) = T(r = T^f) . 



(39) 
(40) 
(41) 



The fourth boundary condition on the gas pressure has to be set at 
the upper boundary of the atmosphere, since we want to calculate 
a consitent value of pgas (i. e., p) at the disc surface. We define a 
constant minimtun value for the mass density 

Pup = p(r = Tup) = const. (42) 

and calculate the corresponding value Pga.s.up at every radial position 
from Tup = T(tup) and pup (see Table 2 for the numerical values of 
Pup and Tup). 

2.3 Numerical solution 

Thanks to the 1 -|- 1-dimensional model, the radial equations decou- 
ple from the vertical structure and can be solved separately. The ap- 
pUcation of the monopole approximation for the disc's self-gravity 
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Figure 1. Temperature dependency of the (a) Rosseland and (b) Planck opacities for certain densities 



requires the enclosed disc mass Mi at radius s to be known for solv- 
ing the radial structure equations (c. f. (6)-(9), (10)). A priori, this 
is only the case at the inner disc radius, where Mj = 0. Due to the 
inner boundary condition, \1' and F tend to zero for s ^ Si (c. f., 
(8), (9)), which causes numerical problems when trying to solve the 
vertical stratification. Hence, we start the calculation close to the in- 
ner boundary, where Mi <^ Mc. First, we solve the radial structure 
equations. With the resulting values of S, ^' and Teff, the vertical 
structure can be calculated numerically in the second step, allowing 
to update the enclosed disc mass (1 1) and to proceed outwards in 
radial direction. 

We use two separate methods to calculate the vertical stratifi- 
cation in the disc and the atmosphere. The disc equations are obvi- 
ously more complicated to solve and as such they are more prone 
to numerical issues like, e. g., steep gradients. We therefore apply 
a Henyey algorithm (Henyey, Forbes & Gould 1964) for solving 
the set of differential equations in the disc. The Henyey method 
looks back on a successful history of applications in stellar struc- 
ture and evolution codes, being able to deal with steep gradients 
by its relaxation method nature. The atmospheric equations, how- 
ever, arc much easier to solve and do not require a powerful, yet 
expensive, algorithm like the Henyey method. We apply a standard 
shooting algorithm to solve the atmospheric structure in a simple 
and quick way. Details about the numerical methods are presented 
in Heinzeller (2008). 

To determine the consistency of the numerical solution for 
the vertical stratification at each radius s, we iterate between the 
Henyey solver for the disc and the shooting method for the atmo- 
sphere. Given an inital guess for p^ff and the boundary conditions 
(27)-(29), the former one provides values for E, h and Teff, once 
the Henyey solver converged to the correct solution. The latter one 
updates the input value pett.cusc from the disc solution by solving 
the atmospheric stratification for the given E, h and Tetf and the 
boundary conditions to peff,atm. The combined solution is accepted 
for 



|Peff,disc — Peff,atm| ^ 6 • min {peff.disc, Peff,atm} , (43) 



with the required accuracy e being defined in Table 2. 



Table 2. Parameters and settings in the numerical model 



Central black hole mass 


Mc 


IM© . . . 


lOOMo 


Accretion rate 


M 


O.OIMe 


..0.15Me 


Standard /3-viscosity parameter 


P 


1 • 10-3 


. . .5 • 10-« 


Corresponding a parameter 


a 


3 • 10-2 


. . . 2 • 10-^ 


Inner disc radius 


Si 


3r-s 




Outer disc radius 


So 


500rs 




Optical depth at upper end of atm. 


'^up 


in-4 




Density at upper end of atmosphere 


Pup 


10-12 




Grid points in s direction 


Ns 


100 




Default grid points in ^-direction 




100 




Maximum grid points in z-direction 


.^z,max 


210 




Max. deviation of disc and atm. sol. 


e 


0.01 





3 RESULTS 

The results presented below were obtained for the set of param- 
eters given in Table 2. Due to the free-fall iimer boundary con- 
dition imposed on (8) and (9), we start the radial calculation at 
s = 2si = 6rs (see also Sect. 2.3). From the values for the sur- 
face density obtained at s = 2si, we can estimate the enclosed disc 
mass for Si ^ s < 2si, finding that its contribution is more than ten 
orders of magnitude smaller than the central mass in all cases. 

We further verified that in all cases the atmosphere contains 
almost no mass, compared to the vertical column of the underlying 
disc. The atmosphere is thin (Zmin ~ h) in in the innermost disc 
region, but expands up to lO/i in the outer regions due to a signifi- 
cantly smaller gravitational attraction towards the disc mid plane. 



3.1 Disc properties of tlie standard disc model 

Our main purpose is to investigate the contribution and efficiency 
of convection in transporting energy and providing viscosity. We 
therefore use a standard setup with Mc = WMq and M = O.IMe 
for which we vary the /3-parameter of the underlying /3-viscosity 



(Me 



1.67 • lOi**g/s(Mc/M0)). We plot the radial structure 



of these discs in Fig. 2 as a function of radius in units of the 
Schwarzschild radius rs = 2.95 • lO^cm. 

Common values for the viscosity parameter /3 are in the range 



of 10" 



. lO"'' (Duschl et al. 2000). To investigate whether the 
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turbulence caused by convection can account partly for the total 
viscosity, we perform disc calculations with /3 = [5 • 10~^; 10~^]. 
We Umit /? to this range for the following two reasons. 

(i) For (i > 10~*, the standard /3-viscosity prescription causes 
the turbulent velocity Uturb./s ~ \/Bsuj to exceed the sound speed Cs. 
In that case, a diffusion limit would have to be introduced (Duschl 
et al. 2000), resulting in an effective decrease of (3 (see Heinzeller 
2008, for a discussion). In the particular example of a lOM© black 
hole accreting at 10% of the Eddington rate, the diffusion limit sets 
infor/3 > 1.4 ■ 10"". 

(ii) For /3 < 10"'\ hardly any solutions can be found for the re- 
quired accuracies and the radial range considered here. The reasons 
therefore will be revealed hereinafter. 

In Fig. 2a, we display the efficiency of convection in the energy 
transport, measured by the dimensionless quantity At each radial 
position, ^ = <^{s, z) is averaged vertically by 




For (3 > 10~*, we find that a significant amount of the total energy 
is transported by convection in the inner part of the disc; close to the 
inner disc radius, C ~ 0.98 for /3 = 10~^. Radiative energy trans- 
port dominates in the outer part of these discs, with a transition zone 
expanding from [10rs;80rs] for /3 = 10"^ to [lOrs; ISOrs] for 
smaller (3. While the curves show a smooth behavior for B ^ 10~®, 
this picture changes when /3 is decreased further. Radial variations 
of C of about 0.1 occur in the case /3 = 5 ■ 10~®, for which we 
also plot a fitting curve. In general, smaller supporting viscosities 
(i. e., smaller values of /3) have little influence on the outer regions, 
while they lead to a significant decrease of the efficiency of con- 
vective energy transport in the inner disc region. 

For a proper explanation of the possible reasons for these vari- 
ations in C wc display further disc quantities in Figs. 2b-g. The 
relative contribution of the radiation pressure to the total pressure 
in the disc's mid plane is shown in Fig. 2b. A comparison with the 
efficiency of convection, described by the quantity nicely con- 
firms theoretical expectations that a strong radiation pressure inside 
the disc drives the convective motion - a simple linear correlation, 
however, cannot be found. We would like to point out that both the 
gas and the radiation pressure do not reflect the instabilities in 

To examine the influence of convective turbulence on the disc 
viscosity, we further display /3avg = VwijvQ, where vq is constant 
for the vertical stratification and 

ph ph 

Vw% = V dz = I {v0 -\- Vcom) dz . 

Jo Jo 

In the low-/3 case, the convective viscosity Ucom becomes three 
times as large as the underlying /9-viscosity. It is important to note 
that although the convective viscosity becomes relatively stronger 
for lower supporting viscosities, its absolute value decreases as 
well. As before, instabilities occur for (3 < 10~^, which are dis- 
played together with the corresponding fitting curve. 

The surface density E increases almost linearly with (3 and 
reflects the variations of ( only very weakly. Since the disc scale 
height h and the density at the disc surface peff are both almost 
unaffected by the value of /3, the increase in E is due to a larger 
internal density in the disc. In all cases, the discs are geometrically 
thin in the outer part, and "slim" in the inner part, with a maximum 
ratio of h/s w 0.2. Contrary to the case of (, no instabiUties are 
found in h for the low-/? case. The density pes at the disc's surface 
shows very similar results for all solutions with a clear decreasing 



trend towards larger radii. A certain irregular structure can be seen 
for all results, an effect of the opacity model, which itself is very 
sensitive to the densities and temperatures in this region of the disc. 
As for the pressure and the disc height, the instabilities in ( are not 
reflected in the density. 

The disc mass Ma increases for decreasing f3, but remains 
completely negligible for all models. We estimate the equality ra- 
dius .Sequ where Mj (s) = Afc by extrapolating the results towards 
larger radii for the [3 = 10"^ disc case. A linear fit to the outer 
region in the log-log plot gives 

Md(s)/Mc = 1 • 10-" (^^j , /3 = 5•10-^ (44) 

which in turn leads to Sequ = 6.7 • 10^ rs. Thus, self-gravity is 
safely negligible in our disc calculations. By means of the radial 
disc equations (6)-(9), this implies the same results for the total 
heat flux F and therefore for the temperature T^tt at the disc surface 
(not shown here), regardless of the value of /3. Furthermore, this 
also implies that the radial variations of C have no effect on the 
temperature profile. In summary, the irregularities of the efficiency 
of convective energy transport are reflected weakly in S, but have 
no influence on the remaining physical quantities. 

3.2 Extended parameter space 

In this section, we extend the disc calculations towards varying ac- 
cretion rates and central masses in order to see how general prop- 
erties and, in particular, the instabilities in ^, depend on the input 
parameters. 

3.2.1 Eddington ratio M / Me 

In the first step, wc investigate the dependence of the results on the 
accretion rate while keeping a constant (3 = 10~^ and a constant 
Mc = IOMq. We perform disc calculations with accretion rates of 
M = [0.01; 0.15]Me. Higher rates are not included, since the discs 
become too thick for the thin-disc approximation to be valid: for 
M = O.ISA/h, the ratio h/s reaches values of 0.3 in the inner disc 
region, while it does not exceed 0.03 for the lower limit M = 0.01 
(see Fig. 3). Furthermore, the same type of radial variations in ^ oc- 
cur for M > O.IMe, which prevent the computations to converge 
for higher values of the accretion rate. For illustration, we display 
them along with the fitting curve in Fig. 3a. For the stable solu- 
tions (M ^ O.IMe), the relative contribution of convection to the 
overall energy transport is smaller for lower accretion rates. This is 
because both the total energy and the angular momentum that have 
to be transported through the disc depend linearly on the accretion 
rate (c. f., (9), (25)). Thus, the standard /3-viscosity is almost large 
enough to accoimt for both requirements when the accretion rate 
is low. Higher accretion rates than O.IMe lead to the same type of 
instabilities of as lower /3-values < 10"^ do for the standard disc 
setup (Mc = lOM©, M = O.IMe). Like in the previous section, 
a higher convective efficiency corresponds to a higher contribution 
of the radiation pressure to the total pressure. 

Figure 3 further demonstrates that the disc mass and the sur- 
face density scale almost linearly with the accretion rate; opacity 
effects modify this scaling law in case of the density pett at the 
disc surface. The contribution of convective turbulence is naturally 
higher the higher the accretion rate is, up to z/com = 3z//3. The con- 
vective zone reaches outwards to 20rs for low accretion rates, and 
to 200rs for high accretion rates, respectively. 
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Figure 3. Solutions for varying accretion rates m = M/M^ with m = Mc/Mq = 10 and /3 = 10 ^ 
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We want to note that the self-gravity of the disc remains neg- 
ligible and therefore the effective temperature scales with Tes oc 
M''^''^'' , as expected from the radial structure equations. As in the 
previous case, the strong variations in ^ are reflected only in the 
viscosity z/avg and the surface density E, though rather weakly. 



3.2.2 Central black hole mass Mc 

Next, we vary the mass of the central black hole in the range of 
Mc = [1; 100]Mq while keeping a constant f3 = 10~* and a con- 
stant Eddington ratio M = O.IMe for the accretion rate. Thus, the 
absolute value of the accretion rate is implicitly scaled with Mc. 
Figure 4 displays selected properties of these disc models. The ef- 
ficiency of convection C, is almost insensitive on a varying central 
mass, with the small differences being to due the disc temperature 
(Figs. 4a,c). Depending on the dissipation rate solely, the effec- 
tive temperature scales with M~°^^ (see (9)). At the same time, 
i^a oc A/c (see (2)). Thus, less energy has to be transported through 
the vertical layers, while at the same time the supporting viscos- 
ity is increased for higher central masses. Radial variations in Q set 
in for the lowest central mass case, Mc = IMq. For completion, 
we would like to add that the ratio Vrng/^fi is almost independent 
of the central mass when the Eddington ratio is kept constant, in 
accordance to the behavior of 

It turns out that both the ratio h/rs and the surface density E 
do not change for varying central masses. For this fact to hold, the 
mass density p has to scale with M^^ , which is reflected nicely in 
Fig. 4b. Then, given that the disc's mass depends only on oc M^, 
the ratio Mi/M^ scales with Mc (Fig. 4d). Note that the calculation 
for the Mc = IOOM0 case terminates at s w 300rs, since the 
density pcff decreases to lO^^^g/cm'', which is the value of the 
upper boundary condition pup in the atmosphere. 

Since the gas pressure (oc Mc~^'^^) decreases more rapidly 
than the radiation pressure (oc M^^) with increasing central mass, 
the ratio prad/ptot is higher the larger the central mass. 



3.2.3 Central black hole mass and Eddington ratio 

Finally, we investigate the case of varying central masses for a 
constant (3 = 10~^ and a constant absolute value M such that 
it equals O.IMe for a lOM© black hole. Thus, the Eddington 
ratio M jMs. scales with M~^. We consider a parameter range 
of Mc = [6.7;100]Mq, corresponding to Eddington ratios of 
M — [0.15; 0.01]AfE. For even lower central masses (i. e., higher 
Eddington ratios), the discs get too thick and also the radial varia- 
tions of ^ become to pronounced to let the calculations converge. 

Figure 5 displays selected properties of these discs, which are 
now controlled by the combined effects of a varying Eddington ra- 
tio and central black hole mass. The efficiency of convection, ex- 
pressed by (" in Fig. 5a, is determined mainly by the Eddington ratio 
and is thus similar to Sect. 3.2.1. In the case of the density at the 
disc surface peff and the temperature T^h (Figs. 5b,c), the two effects 
enforce each other, while the disc mass Md is vastly controlled by 
the central mass and therefore scales as in Sect. 3.2.2 (Fig. 5d). 

The pressure ratio prad/ptot shows an inverse behavior than in 
Sect. 3.2.2, which corresponds to higher convective efficiencies C 
for lower central masses (i. e., higher accretion rates). This inverse 
behavior is due to the fact that, here, the increase in prad oc is 
stronger than the increase in pgas oc p • T for lower central masses. 



3.3 Radial variations in the convection efficiency 

An important point in this discussion is the origin of the instabil- 
ities in for certain disc solutions. We have seen that they occur 
if the underlying viscosity is decreased under a threshold value, 
which itself depends on the parameters central mass and accretion 
rate. Interestingly, these instabilities appear predominantly in and 
and only weakly in the surface density E. The effect on the 
remaining physical quantities is negligible or zero, especially for 
the observables such as the effective temperature. 

Nevertheless, we can understand their occurrence by taking 
a closer look on the vertical structure in the instable zone of the 
disc. We therefore plot the vertical stratification of the ratio of the 
radiation pressure to the total pressure Prad/ptot, the convection ef- 
ficiency C, and the two gradients Vmd, Vad at a radial position close 
to the black hole, s — lOrs (Fig. 6). The data corresponds to the 
case Mc = WMq, M = O.ISA/e and P = W~^, which showed 
significant oscillations of C (Fig. 3). For the horizontal axis, we use 
the heat flux in units of the total flux F, given by the energy 
equation (9). The data is taken from a single solution of the vertical 
structure without any smoothing or averaging. 

For reference, we also display the ratio p,ad/piot at radial po- 
sitions lOOrs and 300rs. Close to the black hole, radiation pres- 
sure dominates over gas pressure, while they equal each other at 
s ^ lOOrs. Further outwards, the disc is gas pressure dominated. 
The average contribution of convection to the energy transport at 
these positions being ^avg = {0.52, 0.1, 0}, we conclude that con- 
vection is radiation pressure driven and by this confirm the results 
of Shakura et al. (1978). In all three cases, the vertical layering of 
the pressure shows smooth curves. Hence, p and T must also adopt 
such a smooth structure and the instabilities in ^ cannot be caused 
by numerical noise in the density or temperature stratification. 

Let us now have a look at the vertical layering of the con- 
vective efficiency ^: we find narrow "convective cells" for small 
Fz ~ 0.1-F-0.5F. These small cells are fluctuating for successive 
iterations (contrary to the extended convective layer between 0.5F 
and F), with the vertical layers close to the mid plane being either 
fully convective or non-convective. The reason for these fluctua- 
tions can be understood from the lower panel of Fig. 6, where we 
display the two gradients Vrad and Vad, which determine whether 
convection takes place in the disc: the quantity C is determined by 
the cubic equation (22), which depends strongly on B, with C — » 
for B and 1 for B ^ 00. The key point is that B reflects 
the Schwarzschild criterium, implying that if the radiative gradient 
is less or equal to the adiabatic one, the stratification is stabilized 
and no convection occurs: 

v„d s; Vad ^ B = 0, c = . 

In the opposite case, even a small positive difference Vrad — Vad 
is multiplied by « 10^° and therefore B > 1 and ( < 1. 
Thus, fluctuations in Vrad — Vad, regardless of being of physical 
or numerical nature, will cause fluctuations in These fluctuations 
can not be seen in the physical quantities, because they occur only 
for small values of Fz and therefore have little effect on the overall 
structure. For a first interpretation of the physics and a relation to 
observable quantities, these instabiUties play only a minor role and 
can be replaced by smoothed values. A further investigation of the 
nature of these fluctuations, however, is necessary in future work. 
For instance, we find that the results showing radial variations in 
also show an inversion in the vertical layering of the density. Con- 
trary to the stellar case, it is not clear whether discs are stable when 
a dominant fraction of the vertical structure is contained within the 
inversion region (Cannizzo & Cameron 1988). 
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Figure 4. Solutions for varying central masses m = Mq/Mq, constant Eddington ratio rh = M/Me = 0.10 and /3 = 10" 



Furthermore, the assumption of hydrostatic equilibrium be- 
comes questionable as soon as the discs are no longer thin and 
the vertical motion is no longer negligible. With a ratio of up to 
h/s ~ 0.3 for small supporting /3-viscosities and for small radii 
s < lOOrs, the resulting discs should be classified as "slim" rather 
than "thin". Interestingly, these cases coincide with the solutions 
showing strong radial variations of C- Thus, in accordance to Jiao 
et al. (2008), waving the assumption of hydrostatic equilibrium and 
including vertical motion might be important for a more detailed 
investigation of these irregularities. 



4 DISCUSSION AND CONCLUSION 

Lower limit on the [5-parameter In the light of the above re- 
sults, we conclude that convection alone cannot account for viscos- 
ity in accretion discs. It requires an underlying viscosity, produced 
by some other process, which is parametrized by V0 in our model. 
The reason for this can be understood from the following line of 
argumentation: convection works towards establishing an adiabatic 
vertical stratification of the disc. Assuming that there exists an ad- 
ditional source of viscosity in the disc, the convective elements are 
decelerated by this inherent friction as well and an equilibrium state 
is established where energy is transported steadily by both radia- 
tion and convection, and where the total viscosity is given by the 
simi of the imderlying and the convective viscosity. If, however, 
the underlying viscosity is too weak, convection is imchecked and 
very efficient in building an adiabatic stratification in the disc with 
Vrad ~ Vad- In such a marginally Schwarzschild-stable state, no 



energy is transported and convection ceases. Thus, convective tur- 
bulence and viscosity vanish. 

In addition, the total viscosity as a result of the vertical in- 
tegration over up + I'conv becomes very small for small values of 
the supporting /3-viscosity. This contradicts the requirements from 
the radial structure equations: the total amount of energy, released 
by the accretion process and given by (9), needs to be transported 
away. Furthermore, viscosity must be present to fulfill the angular 
momentum transport equation (8). Within this argumentation, the 
imderlying viscosity can also be regarded as the "driving force" for 
convection. 

For low supporting /3-viscosities, the density inversion of the 
vertical layering are more pronounced and at the same time, the disc 
is no longer thin. A detailed investigation of the vertical structure is 
needed for a final conclusion about the lower limit for the support- 
ing viscosity. In the limits of this investigation, we conclude that an 
underlying viscosity is necessary and that its minimum value cor- 
responds to a /9-viscosity with [3 ~ 10^ '. This value agrees well 
with laboratory measurements of turbulence induced by differential 
rotation (Richard & Zahn 1999; Richard 2001). 

Influence of central mass and accretion rate Our results show 
that the effects of a varying central mass with fixed absolute accre- 
tion rate are very similar to those of an inversely varying accretion 
rate with fixed central mass. With increasing M /Me, the required 
amoimt of energy and angular momentum that has to be transported 
through the disc increases, leading to larger threshold values for the 
total viscosity. Since convection can only partly account for the re- 
quired increase, the supporting viscosity needs to be larger as well. 
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Figure 5. Solutions for varying central masses m = Mq/Mq with constant M - corresponding to m = M/Me = 0.10 for m = 10 - and /9 = 10" 



Let US consider the case of a varying central mass while the 
Eddington ratio is kept constant. Here, changes in Mc have a strong 
influence on the resulting discs, in particular on the density at the 
disc surface, the effective temperature and the importance of self- 
gravity. Figure 4d suggests that we can assume that the radial scal- 
ing law (44) for Md{s)/Mc holds for higher central black hole 
masses as well. As discussed in Sect. 3.2.2, the ratio Mi{s)/Mc 
also scales with Ale- Thus, the estimated disc mass at s = 500rs 
mcreases from 10""Afc for a stellar mass black hole with lOM© 
to W^'^Mc for a supermassive black hole with 10* Mq, boldly as- 
suming that the extrapolation is valid up to this mass. Correspond- 
ingly, the equality radius Sequ shrinks by a factor 10*. As indicated 
weakly in Fig. 4a, higher central masses in principle allow for lower 
supporting viscosities due to the lower temperatures and a relatively 
stronger supporting /3-viscosity (c. f., Sect. 3.2.2). 

Convective turbulence, differential rotation and magneto- 
rotational instability: a speculative viscosity-mixture Our results 
reveal that disc solutions do only exist if viscosity is also provided 
by effects other than convection. Convection itself can contribute 
significantly to the total viscosity, but needs a driving force to es- 
tablish an equilibrium in energy transport in the vertical direction. 

Here, we parametrize the supporting viscosity by a permanent 
/3-viscosity, where the threshold value of the standard /3-parameter 
depends (weakly) on the central mass and (strongly) on the accre- 
tion rate. For the case of stellar mass black hole accreting at 10% of 
the Eddington rate, we find that « 10~® is sufficiently large, in 
agreement with recent laboratory experiments of rotating Couette- 
Taylor flows. 



In this work, we completely ignore the turbulence created by 
the MRI. Today being regarded as the primary candidate for the 
high viscosity in accretion discs, some aspects still remain to be 
clarified (c. f.. Sect. 1). For example, as detailed in the introduction, 
the question whether the viscosity induced by magnetic effects can 
be translated into an a- or /3-type parametrization is still open. Let 
us assume for the moment that a parametrization is possible. For ex- 
ample, Machida, Nakamura & Matsumoto (2004) investigated the 
case of an accretion disc around a lOM© black hole and found that 
the corresponding a is not constant, but approximately decreases 
linearly with radius: 



a oc exp 



2s/rs 



- 0.99, a 



for s » rs . 



Their results have to be used carefully since the absolute val- 
ues in the fitting formula depend strongly on the disc corona - 
a high-temperature and low-density region, put artificially to pre- 
vent disc material to evaporate (Machida, priv. comm.). Supposing 
that the s^^ behavior of the MRI viscosity is roughly valid, we 
can draw the following picture involving differential rotation, con- 
vection, and magnetic turbulence: in the inner disc region, convec- 
tion and differential rotation with a corresponding /3-parameter of 
~ 10~^ alone do not produce a sufficiently high viscosity for the 
low central mass and/or high accretion rate case. However, close 
to the central black hole, the magnetic turbulence is strong, result- 
ing in a large viscosity due to the magneto-rotational instability. 
In the intermediate disc region, a weaker MRI effect adds to con- 
vection and differential rotation to account for the required total 
viscosity. FinaUy, in the outer disc region, both magnetic effects 
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Figure 6. Vertical disc structure close to the black hole at s = lOrg for 
m = 10, m = 0.15 and /? = 10"^ 



and convection become negligible, but differential rotation is suffi- 
cient in generating the less demanding values of the total viscosity. 
An interesting and important investigation therefore would be to 
combine these three sources of viscosity and to examine whether 
the required viscosity can be generated for a large variety of disc 
parameters. 

In this study, we applied the mixing-length theory to describe 
the convective processes in the accretion disc. Although being ap- 
plied successfully to stellar and accretion disc calculations in the 
past, this theory has several shortcomings like, e. g., the neglect of 
radiative losses and rotation or the imabiUty to derive the anisotropy 



and the mixing-length within the model. We completely ignored 
convection in radial direction, which potentially has significant in- 
fluence on the resulting disc structure through its effects on the ra- 
dial profile of the mass density, for example. In the case of the thin 
discs considered here, however, the radial heat flux is negUgible 
and radial convection therefore not important. 

In the past, alternative theories for convection in accretion 
discs have been proposed, although none of them is fully satisfac- 
tory. For example, Cannizzo & Cameron (1988) investigated the 
importance of convective turbulence in cataclysmic variables and 
compared their results for two different models of convection, the 
mixing-length theory and a self-consistent theory of convection in 
accretion discs (Cabot et al. 1987a,b). Their results showed impor- 
tant differences in the efficiency of convection in generating viscos- 
ity. The interesting project of an investigation of alternative theories 
of convection in our disc model is therefore left as future work. 
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